Linear Continuum Mechanics for Quantum Many-Body Systems 
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We develop the continuum mechanics of quantum many-body systems in the linear response 
regime. The basic variable of the theory is the displacement field, for which we derive a closed 
equation of motion under the assumption that the time-dependent wave function in a locally co- 
moving reference frame can be described as a geometric deformation of the ground-state wave 
function. We show that this equation of motion is exact for systems consisting of a single particle, 
and for all systems at sufficiently high frequency, and that it leads to an excitation spectrum that 
has the correct integrated strength. The theory is illustrated by simple model applications to one- 
and two-electron systems. 



The dynamics of quantum many-particle systems, as 
displayed in electromagnetic transitions, chemical reac- 
tions, ionization and collision processes, poses a ma- 
jor challenge to computational physicists and chemists. 
Whereas the calculation of ground-state properties can 
be tackled by powerful computational methods such as 
the quantum Monte Carlo, [l| the development of simi- 
lar methods for time-dependent properties has been slow. 
One of the most successful methods to date is the time- 
dependent density functional theory (TDDFT), or its 
more recent version - time-dependent current density 
functional theory (TDCDFT).[2] In the common Kohn- 
Sham implementation of this method [3, |3| the formidable 
problem of solving the time-dependent Schrodinger equa- 
tion for the many-body wave function is replaced by the 
much simpler problem of determining N single-particle 
orbitals. However, even this simplified problem is quite 
complex, and furthermore there are features such as 
multi-particle excitations ^] and dispersion forces f6| that 
are very difficult to treat within the conventional approx- 
imation schemes. 

An alternative approach, which actually dates back to 
the early days of quantum mechanics 0, B Q > attempts 
to calculate directly the collective variables of interest - 
density and current. This approach we call "quantum 
continuum mechanics" (QCM), because in analogy with 
classical theories of continuous media it attempts to de- 
scribe the quantum many-body system without explicit 
reference to the individual particles. [lo| 

The possibility of a QCM formulation of the quan- 
tum many-body problem is guaranteed by the very same 
theorems that lay down the foundation of TDDFT and 
TDCDFT.(lll. [l3 | Let us consider a system of particles 
described by the time-dependent Hamiltonian 



energy associated with an external static potential (Vb). 
n(r) is the particle density operator and Vi(r, t) is an ex- 
ternal time-dependent potential. The exact Heisenberg 
equation of motion for the current density operator, aver- 
aged over the quantum state, leads to the Euler equation 



mdtj^{r,t)= - n{r,t)d^[Vo{r) + Vi{r,t)] 

- dyPf,^{V,t). 



(2) 



H{t) = Ho+ / drn{r)Vi{r,t) 



(1) 



where Hq — f + W + Vq is the sum of kinetic energy 
(T), interaction potential energy (W), and the potential 



Here m is the mass of the particles and repeated indices 
are summed over. The key quantity on the right hand 
side of Eq. Q is the stress tensor P^i, (r, i) - a symmetric 
tensor whose divergence yields the force density arising 
from quantum-kinetic and interaction effects. Now the 
Runge-Gross theorem of TDDFT guarantees that the 
stress tensor, like every observable of the system, is a 
functional of the current density and of the initial quan- 
tum state. Thus, Eq. ^ is in principle a closed equation 
of motion for j - the only missing piece being the explicit 
expression for in terms of the current density. 

In recent years much effort has been devoted to 
the theoretic al p roblem of constructing an approximate 

qcm[i1,[i1[i1[I1, 

[itI [isj and several applications have 
appeared in the literature (see Ref. [lil for some represen- 
tative examples). All approximation schemes so far have 
been based on the local density approximation and gen- 
eralizations thereof. In this Letter we derive a new ap- 
proximate expression for P^^(r, t) as a functional of the 
current density for systems that perform small amplitude 
oscillations about the ground-state. The new formula 
is nonlocal, is expressed in terms of calculable ground- 
state properties, and becomes exact in the high-frequency 
limit. 

The Euler equation ^ is conveniently expressed in 
terms of the displacement field u(r, t), which in the linear 
regime is defined by the relation j(r, t) — no(r)9tu(r, t), 
where no(r) is the ground-state density. It is also 
convenient to write the density and the stress tensor 
as the sum of a large ground-state component and a 
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small time-dependent part, i.e., n{r,t) = no(r) +ni(r,t) 
and Pf^„{r,t) = P^^,o(i") + -P/i!/,i(r, Then the time- 
dependent components satisfy the linearized form of the 
Euler equation ^ 

mno{r)dfu = -no(r)VFi(r, t) + .Fi(r, t) , (3) 

where the total force density 

T^A^,t) EE -ni(r,t)apT/o(r) -a^P^^a(r,t), (4) 

is a linear functional of u(r, t). Our approximate expres- 
sion for will be presented in terms of the functional 



E[u] = (iPoHHolM^l 



(5) 



which is the energy of the deformed ground-state j-fAoiu]), 
obtained from the undistorted ground-state {ipo) by dis- 
placing the volume element located at r to a new position 
r -I- u(r, t) . More precisely, we will argue that the force 
density can be represented as 



.F^a(r,t) = - / dr' 



u.(r',t), (6) 



where the second variational derivative of -^[u] , evaluated 
at the ground-state (u = 0) has an exact expression in 
terms of the one- and two-particle density matrices of the 
ground-state. We will show that the representation ^ 
is exact for all one-particle systems and also for many- 
particle systems at sufficiently high frequency. 

Eq. ([6]) can be derived by performinga transforma- 
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a non- 



tion to the "co-moving reference frame" 
inertial frame in which the density is constant and equal 
to the ground state density and the current density is 
zero - and assuming that the wave function in this frame 
is independent of time. This assumption is absolutely 
correct in one-particle systems, where the constancy of 
the density and the vanishing of the current density com- 
pletely determine the wave function. It is also generally 
valid on very short time scales, or for frequencies higher 
than the characteristic energy of single-particle excita- 
tions, because on these time scales it is not possible for 
the particles to "forget" the correlations built into the 
initial ground-state wave function. In all other cases our 
approximation replaces the exact "normal modes" of the 
system by a smaller set of approximate normal modes, 
in such a way that the total spectral weight is conserved. 
We now present a simple derivation of Eq. ([6]), which 
allows us to quickly recognize these facts. 

We start from the linear response of the current density 
to an external vector potential of frequency lu 

j,,{r,u;)= J dr'xp^(r,r',tj)A^,i(r',tj) , (7) 

where (r, uj) is the Fourier component of the current 
at frequency lo and x^j^(r, r', cj) is the current-current 



response function. At high frequency, x^^l' has the well- 
known expansion [20| 

X^u(j,v,uj)^ <5 r - r )5^u + , 8) 

TO m^uo^ 

where the first term (diamagnetic) is frequency- 
independent and 

AV(r,r') = -m\-^a\[[Ho3^.{r^)lWWo) , (9) 

is the first spectral moment of the current-current re- 
sponse function. Now, substituting Eq. ((8]) in Eq. ^ 
and noting that j(r,(jj) = —iujnQ{r)u{r,uj) and that the 
vector potential is related to the scalar potential by the 



equation Ai (r, lu 
der in l/u;'^): 



-, we obtain (to leading or- 



(10) 

This is equivalent to our equation of motion ([3]), with 
•^M.i given by Eq. if and only if 



AV(r,r') = 



(11) 



To show that this is the case we observe that the de- 
formed ground-state is related to the undeformed ground- 
state by the unitary tranformation 



(12) 



Here we have used the fact that the current density op- 
erator j(r) is the generator of a translation of all the 
particles in an infinitesimal volume located at r. Thus, 
the transformation ()12|1 amounts to performing different 
translations by vectors u(r) at different points in space, 
i.e. precisely to deforming the system according to the 
displacement field u(r). Substituting the above expres- 
sion for I^PoM) in the definition of E[u] and expanding 
to second order in u we can easily verify that 

E[u] ~Eo + ^JdrJ dr'u^(r)M^,(r,r')M.(r') , (13) 

which establishes the validity of Eq. (|lip . 

A lengthy calculation allows us to calculate the three 
components of the force density functional arising from 
the kinetic, interaction, and external potential parts of 
the Hamiltonian: = T^]^ + J^^"^ + JC-p°*. The final 

results are 



-j—dud^{nQd^V ■ u) 
1 



+ {2(V^no)u^^ + {dt,no)df,V ■ u 

-f- {dfj,no)d,y\/ ■ u - 25p [{dano)u^a]} , 



(14) 
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I rfr%.(r,r')K(r) 



u,ir')] , (15) 



(16) 



T, 



Here we have introduced the equihbrium stress tensor 

1 



f-LU,0 



1 
2m 



4m 



(17) 

where p*^^-'(r, r') is the one-particle density matrix. The 
interaction kernel K in Eq. (|15p is given by 



/^^.(r,r')=P2(r,r')5^9>(|r-r'|), 



(18) 



is the interaction potential and 
where p*^^' is the two-particle 



where u;(|r — r'|) 

density matrix. Notice that the kinetic force density J-'^^i 
is a semilocal functional of the displacement, involving 
only derivatives up to the fourth order. 

The excitation energies of the system are obtained from 
the solution of Eq. dTU]) after setting Vi = 0. This equa- 
tion defines a hermitian eigenvalue problem with posi- 
tive eigenvalues - the square of the excitation en- 
ergies. The positivity follows from the fact that a de- 
formation of the ground-state wave function must nec- 
essarily increase the energy. The corresponding eigen- 
functions u„ (r) are mutually orthogonal with respect to 
the scalar product (u„,u,„) = / drun{r)um{^)no{'r) — 
if n ^ m. These eigenfunctions must be regarded as 
approximations to the matrix elements of the current 
density operator between the ground-state and the ex- 
cited state in question, i.e. u„(r) ~ i^lno W ' ^^^^^ 

[j]«o(i") = (*n|j(r)|*o)- Even though this is only an 
approximation, it is easy to verify that the sum rule 
J2n^n[j,^]no{r)[ju]on{r') = rR-^ M^^^ir , r') is satisfied by 
the approximate |j]no(r). In this sense our approximation 
preserves the total strength of the spectrum. It is only 
for one-particle systems that the "approximate" [i]no(r) 
becomes exact. Let us now illustrate the theory with two 
simple examples. 

Linear harmonic oscillator. For a harmonic oscillator 
of frequency ujq, external potential Vb(cc) = muj^x^ /2, 

and equilibrium density n^)[x) = — —j^ — , where I = 
(mwo)^^^^, the eigenvalue problem takes the form 



+ [x — 2)u + 3xu — 



2 2 

UJ — UT 



-u = 0, (19) 



where the prime denotes differentiation with respect to 
X. Solving the eigenvalue problem with the boundary 

1 /2 

condition (x)u{x) for oo, we obtain the 

exact excitation spectrum ujn — inujo, where n ~ 1,2, ... 
The corresponding eigenfunctions are Unix) oc i?„_i(x). 
These are indeed proportional to the matrix elements of 
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FIG. 1: (Color online) Unnormalized displacement fields for a 
few low-lying excitations of the two-electron system described 
in the text. The solid line is the ground-state density. Ana- 
lytically we find u„m{x) oc ii'„+m_i[2(a; — a;o/2)] for x ~ xo/2, 
with parity (—1)""^ independent of m. The large value of the 
displacement field for a; ~ does not have a physical signifi- 
cance since the density is exponentially small in that region. 



the current density operator between the ground-state 
and the n-th excited state. 

Two- electron system. Consider a system of two elec- 
trons repelling each other with interaction potential 
1^^^^^ in a one-dimensional parabolic trap of frequency 
ajQ. Due to the separation of center of mass and rel- 
ative variables this model can easily be solved numeri- 
cally, and even analytically in the limit of strong corre- 
lation. We only focus on the strongly correlated limit 
{ujo — *■ 0), since the non-interacting limit turns out to 
be exactly reproduced by our theory. In the strongly 
correlated limit the two electrons become localized near 

1 /3 

X = ±xo/2, where Xq = (2e^/m(^o) is large (see Fig.[T] 
for a plot of the density). The relative coordinate is a 
harmonic oscillator of frequency ujqV^ centered at ixp- 
The center of mass is a harmonic oscillator with fre- 
quency lliq. The exact eigenstates are characterized by 
two non-negative integers n (center of mass) and m (rel- 
ative motion) and are denoted by {n,m). (0,0) is the 
ground-state. The excitation energy associated with the 
state {n,m) is Enm = {n -\- m-\/3)tJo- From the wave 
functions we calculate, without approximations, the dis- 
placement field Unm{x). Some of the results are shown 
in Fig. III). The displacement field of the (1,0) exci- 
tation, which corresponds to a rigid translation of the 
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(n, m) 


exact / 




LO/UJQ 


(1,0) 


1.0 


1.00 


1.00 


(0,1) 


1.732 


1.740 


1.732 


(2,0) 


2.0 


2.643 


2.632 


(0,2) 


3.464 






(1,1) 


2.732 


2.736 


2.732 


(3,0) 


3.0 


3.950 


3.942 


(1,2) 


4.464 






(2,1) 


3.732 


3.965 


3.960 


(0,3) 


5.196 






(4,0) 


4.0 


5.224 


5.217 


(2,2) 


5.464 






(0,4) 


6.928 







TABLE I: Comparison between exact and calculated (appr.) 
excitation energies in the strongly correlated regime. The 
average frequency oj of a group of excitations is calculated 
numerically from the sum rule discussed in the text. Ana- 
lytically one finds Gj'^/ujI = 2 + 3%/3/fc + &k{k - 1)(2 - ^) - 
(— 1)"'(2 — v^)*^, where k = n + m ~ I: these exact values 
are indistinguishable, up to the third decimal digit, from the 
numerical results listed in the leist column. 



Before closing, we point out that the QCM formula- 
tion is appUcable directly to the Kohn-Sham system, in 
which case we do not need the exact ground-state den- 
sity matrices, but only the ground-state Kohn-Sham or- 
bitals and a reasonable approximation for the exchange- 
correlation field. The theory presented here is, in a very 
precise sense, the extension of the well-known collective 
approximation of the homogeneous electron gas to non- 
homogeneous systems and should therefore be useful in 
dealing with collective effects such as multi-particle exci- 
tations and the dipolar fluctuations that are responsible 
for van der Waals attraction. 
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center of mass, is uniform in space, while the displace- 
ment field of the (0, 1) excitation, which corresponds to 
the classical breathing mode, changes sign around the 
origin. The (1,0) and (0,1) modes exhaust the classi- 
cal phonon modes of a system of two localized particles. 
The remaining excitations are fully quantum mechanical. 
Examining Fig. ([T]) one quickly realizes that all the exci- 
tations with a given value of n -I- to and the same parity of 
TO produce the same displacement field, but have differ- 
ent energies. This is a feature of the exact solution that 
cannot be reproduced by any eigenvalue problem with a 
frequency-independent kernel. 

Let us now see what our elastic equation of motion pre- 
dicts for this system. In Table I we present the energies of 
a few low-lying excitations obtained from the numerical 
solution of Eq. (|10p in the strongly correlated regime. We 
see that the energies of excitations such as (1,0), (0, 1) 
and (1, 1), which do not "share" their displacement field 
with other excitations, are very well reproduced by our 
calculation within the accuracy of the numerical work. 
On the other hand, groups of excitations that share the 
same displacement field are replaced by a single exci- 
tation of average frequency, in such a way that the to- 
tal spectral strength of the group is preserved. Indeed, 
it can be proved that the average excitation frequency 
Lj that replaces the frequencies uji of the excitations in 
a given group is given by the sum rule uj'^ = J2i fi^l 
where /; = '^''■i "'(■')-"(■') |.]^g "oscillator strength" 
of the Z-th excitation, u{y) is the normalized solution of 
the eigenvalue problem with eigenvalue uP' , and the sum 
runs over all the excitations in the group. In the last 
column of Table I, we have checked that the sum rule 
is quite well satisfied by the numerical solution of our 
two-electron model. 
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